#ifndef U1_H_
#define U1_H_

#include <cassert>
#include <cfloat>
#include <cmath>
#include <cstdlib>
#include <cstring>
#include <iostream>
#include <vector>
#include <map>

#include "symmetric_matrix.h"

#define NUM_MOVIES 17770
  
/* Timestamp for tracing. */
std::string timestamp();

/* A macro for trace output. */
#define TRACE std::cerr << timestamp() << ": "

/*
 * Functor for comparing with a permutation on a vector.
 * For example, to build perm pi s.t. R[pi[.]] is sorted:
 * vector<int> pi(N);
 * for (int i = 0; i < N; ++i) pi[i] = i;
 * sort(pi.begin(), pi.end(), F_compare_perm<double>(R));
 */
template <typename T> struct F_compare_perm
{
  const std::vector<T>& data;
  F_compare_perm(const std::vector<T>& data) : data(data) { }

  bool operator()(int pi0, int pi1) const {
    return data[pi0] < data[pi1];
  }
};

/*
 * Keeps track of min, max, count, total, mean.
 * Usage:
 * StatCounter trackFoo;
 * ...
 * trackFoo(foo);
 * ...
 * cerr << trackFoo << endl;
 */
class StatCounter
{
  unsigned long long _count;
  long double _total;
  double _min;
  double _max;
  public:
  StatCounter() {
    clear(); 
  }
  void clear() {
    _count = 0;
    _total = 0;
    _min = DBL_MAX;
    _max = -DBL_MAX;
  }

  void operator()(double value) {
    ++_count;
    _total += value;
    if (value < _min) _min = value;
    if (value > _max) _max = value;
  }

  int count() const { return _count; }
  double total() const { return _total; }
  double min() const { return _min; }
  double max() const { return _max; }
  double mean() const { return _total / _count; }
};

std::ostream& operator<<(std::ostream& os, const StatCounter &sc);

/**
 * Sparse histogram.
 * An entry is counted in the k'th bucket iff
 * k*delta <= value < (k+1)*delta
 */
template <typename Value, typename Frequency = unsigned int>
class Histogram
{
  typedef std::map<int, Frequency> Hist;
  public:
  Histogram(Value range, size_t bins)
    : delta(range / bins)
  {
  }

  Histogram(Value delta = 1)
    : delta(delta)
  {
  }

  /**
   * Ordinal of bucket into which value falls.
   */
  int getBucket(Value value) const {
    return (int)floor((double)value / delta);
  }

  /**
   * Puts value into its bucket frequency times.
   */
  void record(Value value, Frequency frequency = 1)
  {
    int bucket = getBucket(value);
    typename Hist::iterator it = hist.find(bucket);
    if (it == hist.end())
      hist[bucket] = frequency;
    else
      it->second += frequency;
  }

  /**
   * Number of measurements in the same bucket as this value, or returns the
   * overflow, if this value is out of range.
   */
  Frequency getCount(Value value) const {
    typename Hist::const_iterator it = hist.find(getBucket(value));
    if (it == hist.end())
      return 0;
    else
      return it->second;
  }

  /**
   * Number of measurements in the map.
   */
  Frequency getCount() const {
    Frequency count = 0;
    for (typename Hist::const_iterator it = hist.begin();
        it != hist.end(); ++it) {
      count += it->second;
    }
    return count;
  }

  /**
   * Number of measurements <= value, where values in the same bucket are
   * indistinguishable.
   */
  Frequency getCountLTE(Value value) const {
    Frequency count = 0;
    for (typename Hist::const_iterator it = hist.begin();
        it != hist.end(); ++it) {
      if (delta*it->first > value) break;
      count += it->second;
    }
    return count;
  }

  /**
   * Percentage of measurements <= value, where values in the same bucket 
   * are indistinguishable.
   */
  double getPercentLTE(Value value) const
  { return (100.0*this->getCountLTE(value)) / this->getCount(); }

  /**
   * Bucket size.
   */
  Value getDelta() const {
    return delta;
  }

  /**
   * Total of recorded data.
   */
  double getTotal() const {
    double sum = 0;
    for (typename Hist::const_iterator it = hist.begin();
        it != hist.end(); ++it) {
      // i * hist[i] can easily overflow an integer, giving the wrong answer
      // when large values are recorded. JLM 20090706
      sum += (it->first)*((double)it->second);
    }
    return delta*sum;
  }

  /**
   * Mean of recorded data, not including the overflow.
   */
  double getMean() const {
    return this->getTotal() / (double)this->getCount();
  }

  /**
   * Standard deviation of recorded data.
   * Here we use the unbiased estimator as a matter of convention.
   */
  double getStdDev() const {
    Frequency samples = getCount();
    if (samples <= 1)
      return 0.0;
    else {
      double mean = getMean();
      double sum = 0.0;
      for (typename Hist::const_iterator it = hist.begin();
          it != hist.end(); ++it) {
        double diff = delta*it->first - mean;
        sum += it->second*diff*diff;
      }
      return sqrt(sum/(samples-1));
    }
  }

  /**
   * Return the smallest value s.t. the fraction of all values less than this
   * value is at least q (e.g. q == 0.9 returns the 90% percentile). Returns
   * 0 if histogram is empty.
   *
   * @param q in [0, 1]
   */
  Value getQuantile(double q) const {
    // Here we exploit the ascending order of keys in the map.
    Frequency samples = Frequency(ceil(q * getCount()));
    Frequency n = 0;
    for (typename Hist::const_iterator it = hist.begin();
        it != hist.end(); ++it) {
      n += it->second;
      if (n >= samples) return delta*it->first;
    }

    return 0; // may reach this if map is empty
  }

  private:
  Hist hist;
  Value delta;

  /**
   * Print in simple one-line format:
   * value frequency value frequency ...
   * with summary stats prepended.
   */
  friend std::ostream& operator<<(std::ostream& os,
      const Histogram<Value, Frequency> &h)
  {
    for (typename Hist::const_iterator it = h.hist.begin();
        it != h.hist.end(); ++it)
      os<<(h.getDelta()*it->first)<<" "<<it->second<<" ";
    return os;
  }
};

/**
 * Natural number sparse histogram.
 */
typedef Histogram<unsigned int, unsigned int> NatHistogram;

/**
 * Minor optimization to avoid call to floor when using natural numbers.
 */
template <> inline int NatHistogram::getBucket(unsigned int value) const {
  return (int)(value / delta);
}

// Arithmetic mean of a vector of numbers.
// Result is undefined if x is empty.
template<typename T, typename S>
T mean(const std::vector<S> &x) {
  T sum = 0;
  for (typename std::vector<S>::const_iterator it = x.begin();
      it != x.end(); ++it)
  {
    sum += *it;
  }
  return sum / x.size();
}

// Variance (N-1 form) when mean of x is already known.
template<typename T, typename S>
T variance(const std::vector<S> &x, T x_bar) {
  if (x.size() <= 1)
    return 0;
  T sum = 0;
  for (typename std::vector<S>::const_iterator it = x.begin();
      it != x.end(); ++it)
  {
    sum += (*it - x_bar)*(*it - x_bar);
  }
  return sum / (x.size() - 1);
}

// Standard deviation of a vector of numbers.
template<typename T, typename S>
T variance(const std::vector<S> &x) {
  if (x.size() <= 1)
    return 0;
  return variance(x, mean<T, S>(x));
}

// Pearson correlation coefficient for two vectors of numbers.
template<typename T, typename S>
T pearson(const std::vector<S> &x, const std::vector<S> &y) {
  assert(x.size() == y.size() && x.size() > 0);
  T x_bar = mean<T, S>(x);
  T y_bar = mean<T, S>(y);
  T num_sum = 0, x_sq_sum = 0, y_sq_sum = 0;
  for (typename std::vector<T>::const_iterator itx = x.begin(), ity = y.begin();
      itx != x.end(); ++itx, ++ity)
  {
    assert(ity != y.end());
    T dx = (*itx - x_bar);
    T dy = (*ity - y_bar);
    num_sum += dx*dy;
    x_sq_sum += dx*dx;
    y_sq_sum += dy*dy;
  }
  if (x_sq_sum == 0 || y_sq_sum == 0)
    return 1;
  else
    return num_sum / sqrt(x_sq_sum*y_sq_sum);
}

/**
 * Add two vectors together component-wise.
 * The vectors must be the same length.
 */
template<typename T>
std::vector<T> &operator+=(std::vector<T> &v, const std::vector<T> &x) {
  typename std::vector<T>::iterator itv = v.begin();
  typename std::vector<T>::const_iterator itx = x.begin();
  for (; itv != v.end(); ++itv, ++itx)
  {
    assert(itx != x.end());
    (*itv) += (*itx);
  }
  return v;
}

/**
 * Scale a vector by a constant.
 */
template<typename T>
std::vector<T> &operator*=(std::vector<T> &v, T x) {
  for (typename std::vector<T>::iterator it = v.begin(); it != v.end(); ++it) {
    (*it) *= x;
  }
  return v;
}

/**
 * Dot product of vectors.
 */
template<typename T>
T operator*(const std::vector<T> &v, const std::vector<T> &x) {
  T sum = (T)0;
  for (typename std::vector<T>::const_iterator it = v.begin(); it != v.end(); ++it) {
    sum += (*it)*(*it);
  }
  return sum;
}

/**
 * L2 norm of a vector.
 */
template<typename T>
T L2_norm(const std::vector<T> &v) {
  return sqrt(v * v);
}

/**
 * Combine multiple keys in a std::multimap using the given combine operator.
 */
template <typename MultiMap, typename F_combine>
void combine(const MultiMap &m, F_combine &combine) {
  if (m.empty()) return;
  typename MultiMap::const_iterator begin = m.begin();
  typename MultiMap::const_iterator end = begin;
  do {
    if (begin->first != end->first) {
      combine(begin, end);
      begin = end;
    }
    ++end;
  } while (end != m.end());
  combine(begin, end);
}

/**
 * Take the average of the values for each key in a multimap and store them
 * in a map.
 */
template <typename Map, typename Iterator>
struct F_multimap_average {
  Map &m;
  F_multimap_average(Map &m) : m(m) { }
  void operator()(const Iterator &begin, const Iterator &end) {
    size_t count = 0;
    typename Map::value_type::second_type sum(0);
    for (Iterator it = begin; it != end; ++it) {
      sum += it->second;
      ++count;
    }
    m.insert(make_pair(begin->first, sum / count));
  }
};

/**
 * Number of distinct keys in a std::multimap.
 */
template <typename MultiMap>
size_t count_keys(const MultiMap &m) {
  if (m.empty())
    return 0;
  typename MultiMap::const_iterator it = m.begin();
  typename MultiMap::key_type k = it->first;
  size_t count = 1;
  for (; it != m.end(); ++it) {
    if (it->first != k) {
      k = it->first;
      ++count;
    }
  }
  return count;
}

bool no_data_for(int movie_id, const symmetric_matrix<float> &S);

#endif // guard

/*
* Copyright (c) 2009 John Lees-Miller
* 
* Permission is hereby granted, free of charge, to any person
* obtaining a copy of this software and associated documentation
* files (the "Software"), to deal in the Software without
* restriction, including without limitation the rights to use,
* copy, modify, merge, publish, distribute, sublicense, and/or sell
* copies of the Software, and to permit persons to whom the
* Software is furnished to do so, subject to the following
* conditions:
* 
* The above copyright notice and this permission notice shall be
* included in all copies or substantial portions of the Software.
* 
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
* HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
* WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
* OTHER DEALINGS IN THE SOFTWARE.
*/

